%% SIMULATION 
% - % step 1: Simulated Rho's from AD 
% -- load the data, save donations AD in rho_AD 
% -- use gendist to generate rho's from probability mass function (relfreq)
% -- define lambdas as the estimated parameters.
% - % step 2: Simulate donations
% -- simulate share that is not affected by deviation cost lambda1, assign
% indicator z=1 to agents that are affected by deviation cost. 
% -- draw deltas from exponential distribution, with mean 1/lambda2 for
% agents with z=1
% -- draw alphas from exponential distribution, with mean 1/lambda3 for
% generous agents
% -- calculate utility levels for all choice options for each agents
% -- each agent picks choice with max utility, save choices in simulated
% donations (done for the three default levels consecutively) 
% - % step 3: save simulated data in CSV file
% - % step 4: Calculate some summary statistics of the simulated data and
% actual data.


clc;
clear all;
format shortE

%% load the matlab data file: prepare data for simulation

load('section4_35.mat');
% use sample with zero contributions added up to 3.5% of sample for 
% illustrative reasons
rho_AD=d.donation(d.default==0);

relfreq= d.relfreq;
%% set seed
rng(4321);

%% take relfreq as underlying distribution: simulate rhos accordingly
N=500000;

P=relfreq'
T = gendist(P, N, 1,'plot');
% rho is T-1 as each element in relfreq corresponds to the value of rho= element -1
rho=T-1;


%% estimated parameters
load('lambdas_hat.mat')

lambda1_hat=x(1);
lambda2_hat=x(2);
lambda3_hat=x(3);
lambda_hat= [lambda1_hat, lambda2_hat, lambda3_hat]
%% Share that acts as if there is are opt out cost (z=0)
n=length(rho);
z=zeros(n,1);
nozero= (rho~=0);
nz=sum(nozero);
% take random draws from uniform distribution between 0-1
% above the value lambda1 --> cost ==1
% random assignment of cost indicator 
% ensures share is about 1-lambda1
z(nozero)=[rand(1,nz)>lambda1_hat];
nd=sum(z);
z=(z==1);

%% draw deltas from the exponential distribution for agents with z==1
delta=zeros(n,1);
mu=1/lambda2_hat;
delta(z) = exprnd(mu, [nd, 1]); %delta exponentially distributed with mean 1/lambda2
disp(mean(delta)); disp(min(delta)); disp(max(delta));
disp(mean(delta(z)));



%% draw the "alpha" cost: the outside option for generous agents
mu3=1/lambda3_hat;
alpha = exprnd(mu3, [n, 1]); %alpha exponentially distributed with mean 1/lambda3
disp(mean(alpha)); disp(min(alpha)); disp(max(alpha));


%% utility levels at different choice options
%V(x=rho)
Vx= (rho.^2)/2;
Ux= Vx- delta;
U0= -alpha;
Ud10= rho.* 10- (10^2)./2;
Ud20=rho.* 20 - (20^2)./2;
Ud50= rho.* 50 - (50^2)./2;
    
%% save data 
% data needed for simulation
save('data_for_simulation.mat', 'rho', 'delta', 'alpha', 'z','Vx', 'Ux', 'U0', 'Ud10', 'Ud20', 'Ud50', 'n');

% data needed for estimation
save('simulated_data.mat', 'relfreq');
%% SIMULATION DEFAULT 10
clc;
clear all;
load('data_for_simulation.mat');


%% individuals' optimal choice

sim_donation10=zeros(n,1);
% if no cost --> always pick rho
sim_donation10(~z)=rho(~z); 

% if endure deviation cost (cost==1)
%  switch to rho if V(rho)-delta>V(d) and Ux> U0
switchtorho=zeros(n,1);
switchtorho(z)= (Ux(z)>=U0(z) & Ux(z)>Ud10(z));
% stay at default if Ud10>=Ux and Ud10>= U0
stay_default=zeros(n,1);
stay_default(z)= (Ud10(z)>=Ux(z) & Ud10(z)>= U0(z));
% opt out completely if U0>Ux and U0>Ud10
optout=zeros(n,1);
optout(z)=(U0(z)>Ux(z) & U0(z)>Ud10(z));

%check the numbers
n_switchers=sum(switchtorho);
n_stayers=sum(stay_default);
n_nocost=sum(~z);
n_optout=sum(optout);
% assign donation to simulated donations
% if switchtorho==1 --> donate rho
sim_donation10(switchtorho==1)=rho(switchtorho==1);
% if stay_default==1 --> donate 10
sim_donation10(stay_default==1)=10;
% if optout==1 --> donate 0
sim_donation10(optout==1)=0;
%% create additional variables (for estimation) and safe simulated data
Delta10=zeros(n,1);
Delta10=((rho.^2)+(10.^2))./2 - rho.* 10;
n=length(sim_donation10);
default10=ones(n,1).* 10;
simd.sim_donation10=sim_donation10;
simd.sim_default10=default10;
simd.sim_Delta10=Delta10

% save final dataset and structure
save('simulated_data.mat', 'simd', '-append');


%% SIMULATION DEFAULT 20
clc;
clear all;
%%
load('data_for_simulation.mat');

%% choices default =20

sim_donation20=zeros(n,1);
sim_donation20(~z)=rho(~z);

switchtorho=zeros(n,1);
switchtorho(z)=(Ux(z)>=U0(z) & Ux(z)>Ud20(z));
stay_default=zeros(n,1);
stay_default(z)=(Ud20(z)>=Ux(z) & Ud20(z)>= U0(z));

optout=zeros(n,1);
optout(z)=(U0(z)>Ux(z) & U0(z)>Ud20(z));


%check
n_switchers=sum(switchtorho);
n_stayers=sum(stay_default);
n_nocost=sum(~z);
n_optout=sum(optout);

sim_donation20(switchtorho==1)=rho(switchtorho==1);
sim_donation20(stay_default==1)=20;
sim_donation20(optout==1)=0;
%%
Delta20=zeros(n,1);
Delta20=((rho.^2)+(20.^2))./2 - rho.* 20;
n=length(sim_donation20);
default20=ones(n,1).* 20;
%%
load('simulated_data');
simd.sim_donation20=sim_donation20;
simd.sim_default20=default20;
simd.sim_Delta20=Delta20
% save final dataset and structure
save('simulated_data.mat', 'simd', '-append');

%% SIMULATION DEFAULT 50
clc;
clear all;
%%
load('data_for_simulation.mat');

%% choices default =50

sim_donation50=zeros(n,1);
sim_donation50(~z)=rho(~z);
switchtorho=zeros(n,1);
switchtorho(z)=(Ux(z)>=U0(z) & Ux(z)>Ud50(z));
stay_default=zeros(n,1);
stay_default(z)=(Ud50(z)>=Ux(z) & Ud50(z)>= U0(z));
optout=zeros(n,1);
optout(z)=(U0(z)>Ux(z) & U0(z)>Ud50(z));

n_switchers=sum(switchtorho);
n_stayers=sum(stay_default);
n_nocost=sum(~z);
n_optout=sum(optout);

sim_donation50(switchtorho==1)=rho(switchtorho==1);
sim_donation50(stay_default==1)=50;
sim_donation50(optout==1)=0;
%%
Delta50=zeros(n,1);
Delta50=((rho.^2)+(50.^2))./2 - rho.* 50;
n=length(sim_donation50);
default50=ones(n,1).* 50;


%%
load('simulated_data');
simd.sim_donation50=sim_donation50;
simd.sim_default50=default50;
simd.sim_Delta50=Delta50;


simd.sim_donation=vertcat(simd.sim_donation10, simd.sim_donation20, simd.sim_donation50);
simd.sim_default=vertcat(simd.sim_default10, simd.sim_default20, simd.sim_default50);
simd.sim_Delta=vertcat(simd.sim_Delta10, simd.sim_Delta20, simd.sim_Delta50);


% save final dataset and structure
save('simulated_data.mat', 'simd', '-append');


%% Write simulated data to csv file
newdata = [simd.sim_donation10, simd.sim_donation20, simd.sim_donation50];

csvwrite('simulated_data.csv',newdata)
%% Summary Statistics Simulated Data versus Original
%{
load('section4_all.mat');

% assign variable names 
default=d.default;
donation=d.donation;

% create subset of observation treatment groups: here excluding observations
% above 300 and zeros
rho_10=d.donation(d.default==10);
rho_10=rho_10(rho_10<=300 & rho_10>0);
rho_10=round(rho_10);
rho_10(rho_10==0)=1;

rho_20=d.donation(d.default==20);
rho_20=rho_20(rho_20<=300 & rho_20>0);
rho_20=round(rho_20);
rho_20(rho_20==0)=1;

rho_50=d.donation(d.default==50);
rho_50=rho_50(rho_50<=300 & rho_50>0);
rho_50=round(rho_50);
rho_50(rho_50==0)=1;


sim_donation10=simd.sim_donation10(simd.sim_donation10>0);
sim_donation20=simd.sim_donation20(simd.sim_donation20>0);
sim_donation50=simd.sim_donation50(simd.sim_donation50>0);

meansim10=mean(sim_donation10);
meansim20=mean(sim_donation20);
meansim50=mean(sim_donation50);
mean10=mean(rho_10);
mean20=mean(rho_20);
mean50=mean(rho_50);

T=table(mean10, meansim10, mean20, meansim20, mean50, meansim50);
display(T);
writetable(T,'mean_data_sim.xls') 
%}
